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ABSTRACT 

We studied a steadily accreting, geometrically thick disk model that selfconsistently takes into account self gravitation of the polytropic 
gas, its interaction with the radiation and the mass accretion rate. The accreting mass is injected inward in the vicinity of the central 
z = plane, where also radiation is assumed to be created. The rest of the disk remains approximately stationary. Only conservation 
laws are employed and the gas-radiation interaction in the bulk of the disk is described in the thin-gas approximation. We demonstrate 
that this scheme is numerically viable and yields a structure of the bulk that is influenced by the radiation and (indirectly) by the 
prescribed mass accretion rate. The obtained disk configurations are typical for environments in Active Galactic Nuclei (AGN), with 
the central mass of the order of 10 7 Mq to 10 s Mq, quasi-Keplerian rotation curves, disk masses ranging from about 10 6 Mq to 10 7 Mq, 
and the luminosity ranging from 10 6 Lq to 10 9 Lq. These luminosities are much lower than the corresponding Eddington limit. 
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CO 1. Introduction 

There is a consensus that a standard model does exist for ge- 
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ometrically thin accretion disks (see e.g. 


Shakura & Sunyaev 


|1973|Lynden-Bell&Pringle 1974 Pring 


e 198 If. Their struc- 



ture is obtained assuming steady accretion and hydrodynamic 
equilibrium. While the viscosity is needed to ensure accretion, 
it can be completely eliminated in favour of the mass accretion 
rate - once the steady disk solution is known to exist - from the 
explicit formulation of structure equations. The luminosity can 
be inferred from the disk geometry and the mass accretion rate, 
but it does not influence the structure of steadily accreting disks. 

In contrast to that, there is no generally accepted model for 
geometrically thick disks. We review just a sample of the re- 
lated literature. PaczyhskT|( |1978[ ) did not solve the full system of 
thick-disk equations, but assumed elements of thin-disk approx- 



imation in determining the vertical structure. Paczyhski & Wiita 
( 1980) imposed some ingredients of thin-disk models onto thick 
disks. The luminosity was obtained analogously to thin disks 
- it was deduced from the already known disk structure. In a 
similar analysis [Abramowicz, Calvani & No bili ( 1980) assumed 
test gas approximation and adapted the mass accretion rate in 
a way suitable to ensure the energy conservation. In a scenario 



discussed injPaczyhski & Abramowicz ( 1982[) the mass accre 



tion rate has been given as part of primary data and the bulk 
of the disk was convective. This description has been applied 
to supercritical accretion, but the selfgravity of the disk was not 
included. |Hashimoto et. al| ( |1993| ) and |Hashimoto et. al| ( |1995| ) 
considered optically and geometrically thick toroidal stars and 
Keplerian disks, with the selfgravity taken into account. They 
observed that the luminosity can be significantly higher than the 
Eddington limit. In these constructions the mass accretion rate 
has been estimated a posteriori. 

The main goal of this investigation is to find a simple self- 
consistent model of steady accretion with radiation and to study 
its structure. We considered a scenario with a geometrically 



thick selfgravitating Newtonian disk (as in Hashimoto et. al 



1995) , in which the mass accretion ra te is prescribed a priori 
(as in Paczyhski & Abramowicz 1982). We assumed, following 



Paczyhski & Abramowicz ( 1982 1, that matter spirals inwards in 
the immediate vicinity of the central disk plane z = 0, where ra- 
diation is generated, by a phenomenological mechanism that is 
inessential for the purpose of this analysis. This radiation inter- 
acts with the gas in the bulk of the disk through Thompson scat- 
tering (i.e., we employed the thin-gas approximation). The struc- 
ture of the bulk is then obtained from the hydrodynamic equa- 
tions (that include a radiation force) supplemented by the energy 
conservation equation. We demonstrate that a self-contained de- 
scription of accretion disks emerges - a simple variant of radi- 
ation hydrodynamics. The whole problem is then analysed nu- 
merically, without any other simplifications. 

Our approach is inspired by the classical accretion model of 
Bondi| ( |1952| ), in which one prescribes the mass accretion rate 
M. The mass accretion rate is not arbitrary, it has to agree with 
other accretion data: the asymptotic mass density p^,, the speed 
of sound flco and the equation of state. In the original model 
of Bondi the selfgravity of gas has been ignored, but an anal- 
ogous construction gives a model that also includes selfgravity 
( |Karkowski et al.| |2P06). As a consequence, the structure of the 
accreting spherical ball of fluid depends on M. It has recently 
been discovered that these models are stable, even if selfgravity 
is taken into account (Mach & Malec 2008 ). 

The content of the remainder of this paper is following. The 
notation is explained and the model is formulated in Sec. 2. 
Sec. 3 deals with radiating test fluids. We explicitly show that the 
equations can be solved by a Born-type approximation scheme, 
which appears to be convergent when the mass transfer rate is 
sufficiently low. Sec. 4. 1 describes an iterative numerical scheme 
that has been specifically invented to solve this problem. Sec. 4.2 
presents the numerical results. The imposed conditions - thin- 
gas approximation, steady accretion and approximately station- 
ary disks - appear quite stringent. We find disk configurations 
with the central mass of the order of 10 7 Mq, Obtained disks are 
geometrically thick, and become even thicker for higher accre- 
tion rates, but the effect of increasing accretion rate and radi- 
ation is not very strong. This is different from features known 
in "Polish donuts" ([A bramowicz, Jaroszyhs ki & Sikora| 
Abramowicz, Calvani & N obili 1980; Paczyhski & Wiita 



1978 



1980 



1 



Patryk Mach and Edward Malec: Accretion and structure of radiating disks 



Paczyriski & Abramowicz , 1982). There is an upper limit for the 
luminosity, and it is significantly lower than the Eddington lumi- 
nosity. These numerical solutions pass the recently formulated 
virial test described in Mach;( [2012| l. 

2. Description of the model 

2.1. Notation and equations 

We assume that a steady accretion disk rotates around a central 
mass M c . The fluid is barotropic, i.e., p = pip), where p is the 
pressure and p is the mass density. Symbols U and <t> will be used 
to denote the velocity of the fluid and the gravitational potential, 
respectively. 

We assume that there exists a phenomenological mechanism 
that produces radiation in the zone occupied by a radial infall 
flow (this coincides with the support of the mass accretion func- 
tion M - see a definition in Sec. 2.2), near the z = plane, 
but we do not specify this mechanism. Pringle remarked that a 
steady disk can be constructed for almost any combination of 
viscosity and radiation process ( Pringle | 1981[ l. We believe that 
one can find - by trial and error - a suitable (point-dependent) 
viscosity near the surface z = that can produce the radiation 
that is obtained in the way described below. 

The radiation is produced with the emissivity per unit mass 
j v and undergoes the Thompson scattering. Our description of 
the radiative transfer follows Padmanabhan ( |2000| l. Neglecting 
the change of the frequency of the scattered photon, we can write 
the scattering cross section as 

i '- =o-<p (fcfcV 

where cr denotes the total Thompson scattering cross section, tp 
describes the angular dependence, and 



j dQtp (k,k') = 1. 



Here k and k' denote appropriate directions in which the radia- 
tion propagates. 

The radiation transport can be described in terms of its in- 
tensity I v . For a stationary process we have 



k-W v 



PJy 
An 



■ CpKly + CpK j doV (£,£') /„(£'). 



(i) 



Here k is the scattering opacity divided by the speed of light c. 
For the fully ionised hydrogen k = o-/(m p c). 
Introducing the radiative flux 



F v = J dQkly (k) 



and integrating Eq. ([T} over the solid angle, we obtain VF V = 
pjy, where we assumed that j v is independent of the direction 
k. The scattering terms that are present in Eq. |l} cancel after 
integration. 

In the following we will use the frequency integrated radia- 
tive flux or the radiation momentum density 



J = I drFy. 
so that 

Vj = p J dvjy 



(2) 



The density and the velocity of the fluid must obey the con- 
tinuity equation 



V (pU) = 0. 



(3) 



The momentum conservation - stationary Euler equations 
with the radiation term - is given by 



(U • V)U = -V<D - -Vp + kj 

P 



(4) 



where <E> is the gravitational potential. The last term describes the 
interaction of gas and radiation in the thin-gas approximation; 
again, only Thomson scattering is taken into account. 

The energy conservation equation states that the radiation 
energy flux j is correlated with the energy flux of the infalling 
gas: 



V ■ (iJp^i + ^U 2 + ojj + V ■ j = 0, 



(5) 



where h is the specific gas enthalpy (dh = dp/p). 

Finally, the gravitational potential is given by the Poisson 
equation 



AO = 4nGp, 



(6) 



where G is the gravitational constant. 

It is now clear that Eqs. ([3]l-(|6]l decouple from Eq. Q and 
any other possible equations describing the radiation transport. 

Let L denote the total luminosity: L — J s dS ■ j, where S is 
the surface of a disk and dS denotes an oriented normal surface 
element. Integration of Eq. Q over the disk volume leads to the 
approximate expression (notice that the neglected term with the 
enthalpy is comparatively small on the boundary) 



(71 



The luminosity depends on the shape of the disk, boundary val- 
ues of the gravitational potential, rotation velocity, and the mass 
accretion rate. All these quantities are intricately related accord- 
ing to the differential equations; the shape itself can be deter- 
mined a posteriori, after solving all of them. 

2.2. Axially symmetric equations 

We will seek axially symmetric solutions of Eqs. ([3])-(|6]l. The 
adopted coordinates are cylindrical variables (r,(f>,z). Here r is 
the cylindrical radius, while R denotes the polar radius, so that 
R = y/r 2 +z 2 . 

The condition of stationarity essentially implies U z = 0; we 
assume that from now on. The velocity vector can be written as 
U = Ud r + cod,),, where we shall demand \ U\ <sc cor. 

Let us define the mass accretion function 

M = -2nUrp. (8) 

Eq. ^ implies that M depends only on the variable z, i.e., M = 
M(z). 

The Euler equations read 

1 i 
<9,-<D + -d r p + Ud r U - co 2 r = Kj r , 

P 



d,<b + -d-p = Kf, 
P 

U (d r (rco) + co) = Kf, 



(9) 
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For barotropes we have Vp/p = V/z, so that V(/i + <t>) = 
— (U • V)U + kj. The expression on the right-hand side has a van- 
ishing curl. Writing this explicitly, one finds that the consistency 
condition is d z (w 2 r) - d z (Ud r U) - -K.(d z f - d r f). We search 
for solutions with small U, and therefore we neglect the term 
d z (Ud r U).lf in addition the rotation velocity u depends only on 
the radius, i.e., co = a>(r), then d z f - d r f = 0, and we can as- 
sume that f = d^, f = d^, where *iF(r, z) is a scalar function. 
It is also convenient to introduce 4* = a:*!*. We shall refer to *¥ (or 
4*) as to the radiation potential. 

Let us take the full divergence of both sides of Eqs. (|9j. After 
trivial rearrangements and neglecting the small term with U one 
obtains 



A^-l-O-T) = l -d r {(D 2 r 2 ). 
Introducing the centrifugal potential 

r , , 2 , 

- I dr r u> (r ), 



(10) 



and integrating equation ( 10 1, we obtain 



h + O + <D C - ¥ = C, (11) 

where C is a constant. Eq. (Hi can be also obtained directly from 
the Euler equations Q. From Eq. (|5]l we have, again neglecting 
terms with U, 



2nr \ r 



(12) 



Combining Eqs. (Ill and ( 12 1, one arrives at 



2nr \ 2 



kM 
2nr 



+ 2ra> 2 



r 2 a>d r (jj\ . 



(13) 



We have to keep in mind that the right-hand side should be eval- 
uated only in the region occupied by the disk. This is a linear 
equation that can be solved iteratively. There exists a unique so- 
lution that is regular in the open space R 3 . 

We observe that the whole problem reduces to Eqs. (jTTJ 
and ( 13 1. Eqs. |6]l and (Hi can be transformed into a nonlinear 
integral equation. The gravitational potential <J can be expressed 
using the Green function formula 



<D(x) 



GM C 
R 



■AnG 



f d\' 
Jv 



G(x - x')p(x'), 



(14) 



where -GM C /R is the contribution due to the central mass, 
V c R 3 is the region occupied by the disk and G denotes 

i.e., 
and 



11 



the Green function of the laplacian in the open space 
G(X) = -1/(4tt|x|). Eq. (fiS) can be combined with Eq. ( 
it is convenient to exploit the fact that the density p is connected 
with the specific enthalpy h by the assumed equation of state. For 
the polytropic equation of state p = Kp T the specific enthalpy is 
h = (KTI(T - 1 ))p T - 1 , and we obtain 



R 



\7lG\ 



T-l 



\ d 3 x' 
Jv 



G(x-x')hr^(x') = C. 



(15) 



In summary, we need the following elements to describe the 
hydrostationary equilibrium of a Newtonian radiating disk: (i) 



the accretion mass function M(z), the rotation curve to(r), the 
central gravitational potential -GM C /R and the equation of state 
p = pip)', (ii) the radiation potential 4*, which can be determined 
from the linear Eq. ( 13 1. It is convenient to write it down using 
the integral equation 



*(*) 



2n Jv r' \2 ) 



(iii) The distribution of the density or the specific enthalpy. This 
can be obtained from Eqs. (|6]l and ( 1 1 1, or from an equation sim- 
ilar to Eq. ( p3] >. 

In all above formulae the rotation law a> = oj(r) is treated 
as known a priori. A popular choice for test fluid solutions is to 
assume a (modified) Keplerian rotation 



GM C 



3/2 ' 



(16) 



where zo is a constant. In this case the centrifugal potential is 
GMjJz- 



Other simple possibilities are the rigid 



rotation a> = const, the rotation law u> = Vn/r, where Vo = const 
is the linear angular velocity, and u> = jo/r 2 , where jo = const is 
the specific angular momentum. We will refer to the last relation 
as to the y'-const rotation law. 

The y'-const rotation is exceptional. In this case we have 

2ru> 2 + r 2 LL>d r co — 0, 

and the radiation equation takes the form 



A^ 



kM 
2nr 



Assuming that 4* is constant at the spatial infinity, we can con- 
clude that 4* = const everywhere. For the proof, notice that it 
suffices to deal with the homogeneous case where 4* — > as 
R — > oo. Multiply the last equation by 4* and integrate over R 3 . 
Integrating by parts and employing the fact that d r M = one 
arrives at 



3 V J-oo V 



The left-hand side is nonpositive, while the right-hand side is 
nonnegative. Therefore 4* = const. That is a torus with a /'-const 
rotation law is not emitting any radiation in r- and z-directions. 
Since the f component of the radiation flux vector vanishes 
identically for the y'-const rotation, we conclude that this rota- 
tion law is not compatible with the radiation. This is consistent 



with the picture emerging from the analysis of Lynden-Bell & 



|Pringle| ( [1974| l; the radiated luminosity balances the energy bud- 
get of the accreting matter and it is accompanied by the shedding 
of the angular momentum. 

The quantity C in Eq. (JTTJ is a free parameter. The bound- 
ary of the disk is defined as a closed two-surface on which the 
specific enthalpy vanishes; thus it cannot be defined a priori; it is 
free. The exception is the test gas approximation without radia- 
tion, where the shape of a boundary is completely dictated by the 
central potential and the rotation curve u>(r). For radiating disks, 
the shape of a boundary also depends on the luminosity, which 
in turn is related to the mass accretion rate. To uniquely define a 
disk, one needs additional information. These questions will be 
discussed in forthcoming sections. 
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2.3. Luminosity and fluxes 

Two of the flux densities are given by f = 3 r v F, / = d^V. The 
third component f can be obtained from the 0-th Euler equation. 

The formula |7]i yields for a disk with the minimal and max- 
imal radial extensions r m and r out , respectively, 



L-J[-.J«/*1»(I 



(tor) 2 + <D 



(17) 



where we employed the condition \U\ <sc ojr. This agrees with 
the formula derived by Lynden-Bell & Pringle (1974) for the 
central star that is co-rotating with the disk, if r m <K r out . The 
local formulae for the flux densities are different. In our case 
the flux density j is defined uniquely, whereas in the standard 
approach it is given up to the total divergence (Lynden-Bell & 
Pringle] 1974} . The accretion flow originates outside of the disk 
and falls onto the central body. It is concentrated close to the 
plane z — 0. The quantity M does not depend on r and thus 
the mass density cannot vanish at the boundary; that means that 
the actual shape of the disk is not well defined near z — 0. 
Nevertheless, in formula ( fTT) we assume that the disk extends 
from a definite exterior cylinder (r out ) to a definite inner cylinder 



3. Test fluid solutions 

Assume that the gas density is low and the gravitational potential 
is dominated by the central term -GM c /7? 3 . It is reasonable to 
expect (and in fact this expectation can be proved) that there are 
solutions of Eq. ( 10 1 that are well approximated by solutions of 
the linear inhomogeneous equation 



Ah = -d r (cj 2 r 2 ). 



(18) 



Consider the rotation law given by Eq. ( 16 1. One can check that 



h = GM C 



R 4 



r 2 + zl) 



(19) 



solves Eq. ( 18 1, with the boundary given by two planes jzj = zo- 
For small zo one recovers the solution of Paczyhski ( 1978| l: 



h = GM C 



Notice that the boundary is not closed. We show below that ra- 
diation causes the closure of the external end of the disk, but it 
is the influence of selfgravity that can make a finite disk. 

3.1. Radiation 

Let the mass of the disk be negligible and the potential be dom- 
inated by the central term <I> = -GM C /R. Assume a modified 
Keplerian rotation curve ( [To} . The solution of Eqs. ( fT0| ) and ( fT3] l 
can be found perturbatively, treating the mass accretion term as 
a perturbation. 



The zeroth-order term ho is given by Eq. ( 19 1, that is 



GM C 



<-+z, 



0) 



(20) 



and the hi term is 



■X- 



//,(x) - | d\iG(xi-x)^-d ri ^> + h + ^oj 2 r 2 j 

kM 



( d 3 xiG(xi-x)^^GM c x 
Jv 



1 



Higher order terms h k {k -2,3,. . .) are defined by 



(21) 



h k (x) 



-I 



kM 



d x k G (x k - x) - — d n h k -\ (x k ) . 
v 2nr* 



(22) 



Here, for £ = 1,2, ...,x k = {r k ,z k ,<p k ); r k , z k and <f> k are the 
cylindrical components of x k , and the integration volume ele- 
ment d 3 x k equals d 3 x k = dr k dz k d<p k r k . The unlabelled quantities 
r, z and <p are reserved herein and in what follows for the com- 
ponents of the vector x = (r, z, 4>). 

The specific enthalpy corresponding to the full solution of 
Eqs. ( 10 1 and ( 13 1 is h = ho + h\ + 6h, where 5h = hi + hi + . . . 
It is easy to check that all individual functions hi, hi, /13, ... are 
non-positive. Assume furthermore that the disk extends from the 
inner cylinder of a radius r m to an outer cylinder of a radius r out . 
Then it is a simple exercise to show that 



\d n h k \ < 



\h k 



(23) 



for any k — 1,2,... 

Stationary disks can exist when the above expansion scheme 
is convergent. The application of standard criteria for conver- 
gence of series leads to a bound onto the mass accretion rate M 
and on the total luminosity. These statements can be proved quite 
generally. We would like to point out that the rotation law ( 16 1 is 
not compatible with the geometry of a closed disk. While it al- 
lows for a closure at a finite distance from the centre, it should be 
modified in the interior region to give a disk configuration. Let 
us point out, however, that the inclusion of self-gravity should 
allow for a closed disk configuration. 

3.2. Radiation and the boundary of the disk 

The boundary of a disk is a smallest surface where the specific 
enthalpy vanishes: h = ho + h\ + 5h = 0. In the absence of 
accretion, i.e., when M — 0, we have h k = for k = 1,2, .. . 
Thus at z = Zo there exists a horizontal part of the boundary. 

If the mass accretion term M is positive, then it is clear from 
the inspection of Eqs. (21 1 and (22i that h\ and 5h are non- 



negative. This means that the boundary is pushed inward to a 
location with a value \z\ < zo- 

The directional derivative of the specific enthalpy h vanishes 
along the boundary: dh — 0, dh — d r hdr + d,hdz. Notice that 
Kj r = d r {h\ + 6h) and kj, = d z (h\ + 6h). This leads to the equa- 
tion 



dz 
dr 



Kjr + 9 r ho Kj r + d r h 



Kj z + d z ho 



(24) 



Values of j r and j z are positive on the upper face of the disk 
(z > 0) and negative on the lower face (z < 0). For vanishing 
radiation, Eq. ( |24| > reduces to the well-known formula dz/dr = 
(u> 2 r - d r <I>) /a<D. 
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3.3. A family of analytic solutions 

Assume M = F(r)5(z), where F(f) = F — const > for 
r m < r < r out and F(r) — for r < r m and r > r out . The 
delta-like distribution can be easily handled analytically, and a 
close inspection shows that some aspects of the disk solutions 
weakly depend on the specific form of the mass accretion func- 
tion M. We assume that F « J dz(x>r 2 p, because the radial drift 
should be negligible in comparison to the rotational motion. For 
r m < b < r < r out the rotation curve is assumed to be given 
by Eq. 



Eq. (20 1. With this assumption one can obtain the cusp-like end 
of the external part of the disk. This rotation law does not allow 
to close the internal part of the disk. To achieve this, we have 
to modify the rotation curve in a transient region just above r m . 
This transient zone is not particularly important, since it does not 
seriously impact the structure of the disk nor its luminosity. For 
instance, one can take in the interval r m < r < b, 



GM C 



3/2 



2 + or- 1 + 



(25) 



The constant F describes the intensity of the radial baryonic 
drift, and it can be found from the condition that h(r out ) = 0. 
The parameter a can be specified by the demand that h(ri n ) = 0. 
Notice that the rotation curve ( 25 1 yields the specific enthalpy 



h = GM C 



1 2 - (r/b) a 
R 



4 



r2 + 4 ) 



in the region {r m , b). It is continuous everywhere, albeit it is not 
differentiable air - b. 

Let us define the elliptic function 



E(f,r k ,z) 



X 

4^ Jo 



d(p k G (x - Xjfc) 



^jr 2 + r 2 k +z 2 - 2r k r cos (f> k 



(here k — 1,2,...). The term hi is now given by 
kFGM q 



h\{r,z) = -- 



An 



f 



dr\r\E{r,r\,t) x 



3zl 



(26) 



2/r 



X 



kfgm c r h , 

dr\E(r,r\,z) X 



xd r 



(i- (! + «)(*)" 4(1-0)1 



One can specify the parameters b and a of the transition so that 



\3/2 



the second integral in Eq. ( 26 1 can be neglected, at least for thin 
disks. 

The next expansion terms h^s (k = 2, 3, . . .) are given by 



h k (x) 



kF 
~2~n 



£ 



dr k E(r,r k ,z)d rt h k _i 



We can use estimate ( |23) to obtain the inequality 

£ dr k E(r,r k ,z) J sup \h k -i\ 



kF 

sup |^| < sup 



2nn, 



(27) 



The elliptic function E(r,r k ,z) is integrable on any finite disk. 
The inequality ( p7| means that the series expansion is convergent 
for sufficiently low values of the product 



kF 



16 i, which yields the h^ enthalpy term according to 2nr n 



£ 



dr k E(r,r k ,z). 



It suffices that the mass current constant F is small enough, 
which means - notice the conservation law ( 17 1 - that stationary 



disks cannot have arbitrarily high luminosity. 

The mass accretion rate influences the disk geometry, as seen 
from the following discussion. At the boundary we have h^+h\ « 
0; thus (keeping only the leading term of hi) we have 



1 



1 



R 4 



r^+zi 



kF 
4^ 



dr\r\E(r,r\,z) x 



3d 



for r > b and 
1 2 - (r/b) a 



R 4 



r 2 +z\ 



kF 
An 



(-7 + 4)' i 

X 



(28) 



dr\TiE{r,ri,z) x 



3z,1 



5/2 



(29) 



for r < b. Assuming that the outer part of the disk extends up to 
r out , one obtains the value of F 



kF 



4tt 2r mt I' 



(30) 



where 



dxxE ( 1 , x, 0) 



2/^2 
out 



5/2 



One can check that 

dr k E (r out , r k , 0) > y(r in /r out ) 



£ 



dr k E(r,r k ,z) (31) 



for small zo- Here y(ri n /r out ) is a coefficient, with values de- 
pending on the ratio r; n /r out , ranging between 0.01 and 1 for 
In/rout = 10~ 9 , . . . , 1. It is found numerically as the largest co- 
efficient ensuring that the inequality 



X' 



dxE(l,x,Q) > y(r in /r out ) 



£ 



dxE(r/r 0Ut , x, 0) 



is satisfied. Inequality ( 3 1 1 combined with Eq. ( 30 1 yields 



kF 
2nn, 



£ 



dr k E(r,r k ,z) < 



y(rm/r ut)rinr 



Obviously, for thin disks the right-hand side is much smaller than 
one; this implies the convergence of our approximation scheme. 



5 



Patryk Mach and Edward Malec: Accretion and structure of radiating disks 



The consistency condition F <K J dzu>r 2 p, discussed earlier, can 
be checked only a posteriori, after solving all equations. 

Eq. (29 1 allows, putting r = r; n , to specify the parameter a 
that appears in the rotation curve in the transient zone (r; n , b). 
By a proper choice of b, r m , r out and zo one can always achieve 
a <K 1, which gives the total luminosity close to the value cor- 
responding to the Keplerian rotation curve. Formula (17i now 
yields 



FGM r I 1 



1 

r out 



Other quantities, including the gas density, can be found from 
the Euler equations and the radiation transport equations. 



4. Heavy selfgravitating disks 

4. 1 . Numerical methods 

Solutions corresponding to heavy selfgravitating disks can be 
obtained numerically by solving Eqs. (|6|, ( [TT] ) and ( 13 I. In the 
following we assume the polytropic equation of state p = Kp r , 
so that h = {KT/(T - l))p r ~'. The rotation law is fixed up to 
a multiplicative constant. Numerical examples will be given for 
Keplerian-type rotation u ~ r~ 3 / 2 . 

Assume that the disk spreads up to R = R out at the equa- 
torial plane, and that the maximal density of the gas is p mdx . 
The quantity u = GR 2 ut p m - dx has the dimension of the potentials 
<E>, <t> c and 4*. It can be used to define following dimensionless 
quantities O = O/m, O c = <5 c /m, 4* = "¥/u, K = Kp^/u, 
ib = toRoutl V"- I n similar fashion we introduce p = p/p mdx and 
M c - M c /(p m . dx Rl ut ). Dimensionless spatial coordinates are de- 
fined as x = x/R out . The relevant equations of the model can be 
rewritten in the form 



Kr/(r - l)p 1_1 + <5 + <5 C - Y = C, 
AO = 4-np, 
AT = S, 



(32) 
(33) 
(34) 



where C = C/u, A is the laplacian with respect to the rescaled 
coordinates x, and 



kM 
2n? 



-(d,^ + 2~ru> 2 + r 2r 
? v 



The system of Eqs. (32i and (|33J) with *P = 0, i.e., with 



out radiation, corresponds to a simple model of a rotating star 
or a toroid. There is a long tradition in the numerical analysis 



of these equations (see e.g. Stoeckly 



1965 Ostriker & Mark 



196^|EriguchI|[T978l[Enguchi & Muller||1985| >. In particular, it 
is known that they can be solved iteratively starting with some 
initial guess of the density distribution. The very fruitful iter- 
ative scheme of Ostriker & Mark (1968) based on the Green 
function expression for the gravitational potential is known as 
the self-consistent field method. It has been used successfully 
by many authors and in many variants, including computation 
of general-relativistic solutions (cf. Clement| 1974| Blinnikov 



1975 Komatsu, Eriguchi & Hachisu 1989 Nishida, Eriguchi & 
Lanza||1992[). An anal ytic expansion scheme has been also pro- 
posed by Odrzywolek (2003). Since in our case Eq. (34i for the 
radiation potential is also Poisson-like, the numerical method of 
this paper follows the self-consitent field pattern. 



Knowing the density distribution, one can compute the grav- 
itational potential <E> from a rescaled version of Eq. ( 14 1, i.e., 



O(x) 



M c „ 

~ +4tt 
R 



d 3 xG(x - x)p(x). 



(35) 



Eq. p5\ is clearly not suited for a direct numerical implementa- 
tion due to the singularity at x = x present in the Green function 
G{x - x) = — 1/(4tt|x - x|). This difficulty can be avoided by 
a regularisation procedure. The standard trick is to expand the 
Green function in terms of spherical harmonics. Assuming ax- 
ial symmetry, and spherical coordinates (R, p = cos 6, (j>) we can 
write 



, p) = - ° - An V P 2j (ft) dRR 2 F 2j (R, R) x 

f 

Jo 



where 



Fj(R,R) = 



x dpP 2j (p)p(R,p), 



(36) 



and Pj denotes the j-th Legendre polynomial. From the numeri- 
cal point of view the regularisation of the discretised integral in 
expression ( 35 1 consists of taking a finite number of terms in the 
series in Eq. (36i. 

The equation for the radiation potential T* can be solved in a 
similar fashion. Assuming that § is known, we can obtain < P as 



°° /-»oo 

y(R,p) = -Y.P20) dRR 2 F 2j (R,R)x 



x dpP 2j (p)S(R,p). 



(37) 



Knowing *F is necessary to compute S , which in turn is needed 
to find S . We demonstrate that in the case investigated here an 
iterative process can be applied that eventually yields an appro- 
priate solution both for § and T*. The expression for S is valid in 
the interior of the disk, i.e., where p + 0. The integral in Eq. ( [37} 
has to be evaluated only over this region, hence p is implicitly 
assumed to be known. 

Conversely, once 5> and T* are known, we can compute the 
density p directly from Eq. |32i. This equation, as it is written, 
allows negative values of enthalpy. In this case one has to mod- 
ify the physical boundary of the disk. We will only search for 
p in domains where the enthalpy given by Eq. ( |32] > is positive; 
otherwise we will assume p — 0. 



The strategy of solving Eqs. (32 1, (33 1 and (34i is as fol- 
lows. We fix the value of the central mass M c , and the ratio of 
the inner and outer radii of the boundary of the disk Ri„/R out . 
In the next step we assume temporarily that T* = 0, and con- 
struct an initial density distribution in the form of a toroid with 
assumed ratio R m /R ollt and the maximum density p = 1, For this 



toroid, the gravitational potential can be obtained form Eq. ( 36 



and the new density distribution can be computed from Eq. ( 32 
This procedure is iterated until a desired convergence is reached 
After establishing a hydrostatic structure of the selfgravitating 
disk, we "switch on" the radiation. We fix the function kM{z) and 
obtain the radiation potential coming from the already found hy- 
drostatic structure. This can be performed by solving iteratively 
Eq. ((37]). Finally, all three Eqs. ((36]), ([37]), and Q are solved 
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Fig. 1. Density p obtained for solution (a). The plot shows a sec- 
tion of the upper hemisphere in a meridian plane. The density is 
greyscale-coded. 

iteratively in the same fashion. This eventually produces a con- 
figuration with the hydrostationary structure that is influenced 
by radiation and mass accretion. 

A subtle point in this procedure is that the appropriate val- 
ues of constants C and K, as well as a multiplicative constant 
appearing in the rotation law, are not known a priori. One has to 
perform a kind of renormalisation of these constants during the 
iterative process described above. From the condition that the 
inner and outer edges of the disk are located at R[ n and R out , re- 
spectively (i.e., the density must vanish there), we can compute 
the values of C and the multiplicative constant in the rotation 
law, provided that the potentials <1> and T are known. Expressing 
the rotation law as Cj{?) = a>f(f) with u> denoting the appropriate 
multiplicative constant, we can write 

— (5(r= 1,2 = 0)-¥(r = 1,2 = 0) 

-6(r = RJRouul = 0) + <F(r = R m /R out ,~z = 0)) . (38) 



oj 2 = 



The value of C can be now obtained from Eq. ( |32[ ) taken at /?;„ 
or R out with p — 0. The constant K is given a value that en- 
sures that the density corresponding to the found maximal value 
of enthalpy h is precisely p max . These computations have to be 
repeated each time a new density distribution is computed. 

Convergence properties of our numerical scheme are depen- 
dent on the spatial resolution of the numerical grid. In this paper 
we follow the optimisation techniques for the computation of 



integrals appearing in Eqs. ( 36 1 and ( 37 1 that are described by 



Miiller & Steinmetz ( 1995 1. They allow one to achieve spatial 
resolutions of the order of 5000 x 5000 on a parallel computer 
consisting of 64 processors, with a total computing time of sev- 
eral minutes only. We truncate the expansion in Legendre poly- 
nomials Pi in Eqs. ( 36 1 and ( 37 1 around Z max = 20. The conver- 



gence is mainly controlled by computing the maximum norm of 
the two subsequent density distributions in the iteration scheme, 
i.e., we compute e = ||p' +I -p* , where the index i numbers 
subsequent iterations. We continue to iterate until e reaches the 
level of 10~ 6 , or even 10 8 . Convergence in 3> and "T is controlled 
in a similar fashion, the difference being that the maximum ab- 
solute values of these quantities are not known a priori. 

4.2. Numerical examples 

We discuss two sample solutions (a) and (b), obtained with 
the method described in the preceding section. They are com- 
puted for the poly tropic exponent F = 5/3, and the central 
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Fig. 2. Part of the gravitational potential due to the accretion 
disk only, i.e., <t> + M c /R, obtained for solution (a). The plot 
shows a section of the upper hemisphere in a meridian plane. 




Fig. 3. Radiation potential obtained for solution (a). The quantity 
log 10 |*F| is plotted in greyscale. The graph shows a section of the 
upper hemisphere in a meridian plane. 



mass M c = 2. We assume the quasi-Keplerian rotation law 
a> = (l>/r 3 l 2 . We refer to this rotation curve as "quasi-Keplerian" 

because u> + yjM c > although the actual difference between a> 

and V^c appears small. Notice, for the subsequent discussion, 

that setting a> = V^c is equivalent to co(r) = yGM c /r 3 '' 2 . We 
take kM distribution in the form 



kM = Aexp(-(50z) 2 ). 



(39) 



Solution (a) is obtained for Ri n /R on i = 10 4 and A — 10 4 . 
Solution (b) is computed assuming R m IR ut = 10~ 3 and A = 
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10~ 2 . The density p, gravitational potential of the disk <b + M c /R, 
and the radiation potential *F, obtained for solution (a) are shown 
in Figs. [T| [2] and [3] respectively. We do not display solution (b), 
because its shape is not much different from that of solution (a). 
The obtained value of a> is almost that of a strictly Keplerian 

rotation law. Expressing a> as cb — (1 + a) ~<J~M^ we obtain a * 
4.7 ■ 1(T 6 for solution (a) and a » 4.0 • 1(T 5 for solution (b). 
This can be intuitively understood by the careful inspection of 
Eq. ( 38 1; for elongated disks (i.e., for low values of Rm/Rout)> the 
total gravitational potential is dominated by the divergent term 
-GM C /R near the centre. 

Other physical parameters of the solution depend on the 
choice of R out and p max , as well as on the actual values of 
constants such as k and G. Possible values of R oal and p max 
are restricted by the assumptions that were made when deriv- 
ing equations of the model. The assumption of the thin-gas 
approximation requires that the mean free path of the photon 
A - m v l(crp) = (cKpY 1 should be comparable with the size of 
the entire configuration. The radial velocity \U\ should be lower 
than the angular one ra>. Finally, the temperature of the disk 
should be high enough so that the assumption that the gas is 
fully ionised is justified. The radial velocity U can be computed 
from Eq. ([H} based on the density distribution and the assumed 
mass accretion rate M. Assuming the perfect gas approximation 
and a value of the mean molecular weight of the gas p., we can 
also find the distribution of the gas temperature 

T = ppm p /{pk B ) = Kpm p p r ~ l /fc B , 

where k B denotes the Boltzmann constant. For the fully ionised 
hydrogen one can take p = 1/2. The maximal temperature in the 
disk can be computed as 

r ma x = upm p K/k B . 

We take R out = 10 parsecs for both solutions, and as- 
sume p max = 10~ 18 g • ctrT 3 and p max = 10~ 17 g • cirT 3 for so- 
lutions (a) and (b), respectively. This gives the central mass 
M c = 2.9544 • 10 7 M© for solution (a). The central mass for 
solution (b) is higher by a factor of 10. Thus, if solutions (a) and 
(b) were to be applied in an astrophysical context, they would 
serve as models of accretion disks around ultramassive galactic 
black holes. 

Other physical quantities characterising solutions (a) and (b) 
are as follows: For solution (a): M fluid = (3.18720 ± 0.00049) • 
10 6 M o , L = (3.1101 ± 0.0092) • 1O 6 L , r max = (3.73508 + 
0.00031)- 10 4 K. For solution (b): M fluid = (3.20779 + 0.00056)- 
1O 7 M , L = (2.84443 ± 0.00088) • 10 9 L o , r max = (3.72160 + 
0.00036) • 10 5 K. The error estimates given here were obtained by 
comparing solutions computed on numerical grids of different 
resolution (spanning the range of 2000 x 2000 to 4000 x 4000), 
with different numbers of the Legendre polynomials (/ max = 16 
to / max = 22), and different convergence level (e = 10~ 6 up to 
e = 10~ 8 ). 

The estimate of the mean free path of a photon A - evaluated 
for the maximal gas density - is of the order of 1 parsec for 
solution (a), which is roughly the maximal vertical height of the 
disk. For solution (b) this estimate is 10 times smaller, but it 
is clear from the shape of the disk displayed on Fig. [T]that the 
optical thicknes is smaller than one, i.e., 



t = sup 

r£[R m ,R a 



r 

1 Jo 



crp(r, z) 



dz<\. 



Numerical correctness of our solutions can be tested using 
the following virial theorem. Define £p 0t = J d 3 xp(<t> + <E>Kep)/2, 
Ekin = J d 3 xp\U\ 2 /2, E theim = J d 3 x3p/2, and E = k J d 3 xpx ■ j, 
where <t>K ep = -GM C /R is the Keplerian potential of the point 
mass only (note that <£> Kep enters the formula for E pot twice - <D 
is the total gravitational potential, i.e., the sum of <J>K ep and the 
contribution from the accretion disk). It can be shown (Mach 
|2012[ ) that for a stationary configuration consisting of a point 
mass M c and the fluid satisfying the Euler equation Q the fol- 
lowing virial identity holds: 

£p 0t + 2£ , ki n + 2E t h erm + E = 0. 

It differs from a standard formulation of the virial theorem by 
the presence of the point mass term in £ pot , and the radiation 
coupling term E. The virial check of a numerical solution can be 
performed by computing 



6 V = |£pot + 2£'kin + 2Ethenn + E\/\E pot \. 



(40) 



We have obtained e v w 10~ 8 for all numerical solutions. All 
constituent terms appearing in ( |40| were always much greater 
than this value. The following quantities are roughly the same 
for both solutions (a) and (b): £kin/|£potl ~ 0.478, E tam /\E pot \ « 
2.22 • 10- 2 , (J d 3 xp®/2) /£ pot * 0.525, (/ <i 3 xp(I>Kep/2) /E pot * 
0.475. The ratio £/|£potl equals approximately 3.08 ■ 10~ 6 for 
solution (a) and 2.81 ■ 10~ 4 for solution (b). That confirms the 
validity of our results. Clearly, our solutions pass the virial test 
surprisingly well. That is not very common, but not exceptional 
either (see e.g. Axenov & Blinnikov 1994 1. 



5. Summary 

We have formulated a consistent model of a selfgravitating disk 
with steadily accreting matter. The radiation emitted by the disk 
interacts with the gas by Thompson scattering (thin-gas approx- 
imation). The most interesting feature of our solutions is that the 
accretion mass rate flux density is concentrated in the equato- 
rial plane z — 0. The slow radial drift of gas is observed in the 
equatorial plane and the bulk of gas remains approximately sta- 
tionary. It appears that the conservation laws of the energy and 
the momentum together with the assumption of approximate sta- 
tionarity suffice to obtain the structure of the disk. The emissivity 
index of accreting matter can be deduced after solving equations 
of the model. The approximation of stationarity demands that the 
radial inflow speed of matter has to be negligible compared with 
the rotational velocity of the gas in the disk. The secular change 
in the mass of the central accreting object should also be negli- 
gible. These two assumptions have been verified post factum for 
the numerical solutions presented in this paper. 

The mathematical description of the radiating disk reduces to 
a pair of elliptic partial differential equations. They can be solved 
iteratively in a way similar to that routinely used in the litera- 
ture when finding selfgravitating equilibria of non-radiating gas 



(cf. Ostriker & Mark 1968; Eriguchi & Miiller 1985 Nishida, 



Here zo = Zo(r) is the disk height at the cylinder of radius r 
UMihalas & Weibei-Mffialasl[l98% 



Eriguchi & Lanza 1992| l. In our case each iteration step con 
sisted of finding new distributions of the density, the gravita 
tional potential and the radiation potential. One of the main re 
suits of this paper is that this procedure numerically converges. 
Analytic results can be obtained in simplified cases. We investi- 
gated the influence of the emitted radiation onto the disk struc- 
ture in the test fluid approximation. The interesting conclusion 
is that approximately stationary solutions do exist only when the 
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mass accretion (and thus the luminosity) is not too large. This in- 
tuitively well-understood feature of solutions has been revealed 
also in our numerical analysis of heavy selfgravitating disks. 

We assumed thin-gas approximation, slow radial flow of 
matter, and approximate stationarity. These conditions appear 
quite restrictive - in the parameter space of solutions there is 
only an island of parameters for which the above assumptions 
can be satisfied. They lead to solutions with essential features 
that agree with the conditions encountered in some AGNs. We 
found solutions with the central mass of the order of 10 7 Mq to 
10 8 M© and disk masses of the order of 10 6 M© to 10 7 M©. The 
luminosity varies between 10 6 L© and 10 9 L©. We assumed the 
quasi-Keplerian rotation law u> ~ r~ 3/2 , but the rotation appears 
in fact to be Keplerian in the examples described in this paper. 

It is interesting that the luminosity in numerical examples 
has been much lower than the Eddington limit. This differs from 
the results of [Paczyhski & Abramowicz (1982) or Hashimoto] 
et. al ( |1995 1, where supercritical luminosity has been discovered. 
This discrepancy has a physical explanation - we describe disks 
in the thin gas approximation, while the quoted works dealt with 
the disks in thermal equilibrium or with convection transport. 

A future investigation of our model can be aimed in two di- 
rections: the study of stability of disks and the formulation of a 
general-relativistic version. Unlike for the standard models, the 
structure of our Newtonian disks is dependent on the accretion 
flux and radiation. We have discovered that Bondi-type, spher- 
ically symmetric solutions are stable also in the selfgravitating 
regime (Mach & Malec 2008). The Bondi accretion models are 
spherically symmetric in contrast to accreting disks, but they 
share with our model the property that their structure depends 
on the accretion, and that can have a stabilising effect also in the 
nonspherical case. Thus we expect that the solutions discussed 
in this paper are stable. 

There are two extreme classes of general-relativistic radiat- 
ing accretion disks. Disks characterised by the size of an inner 
boundary exceeding 2GM c /c 2 (quasi-Newtonian case), where 
M c is the central mass, should have a structure similar to those 
presented in this paper. However, their stability properties can 
still be different owing to the influence of gravitational radia- 
tion. Disks that are inherently relativistic, with the inner bound- 
ary located close to the minimum stable surface (6GM c /c 2 for 
the Schwarzschild solution), can be significantly brighter than 
their Newtonian counterparts - again, this expectation is based 



on our experience with Bondi-type solutions ( |Karkowski et al. 
[20061 palec & Rembiasz1|20T0) i. 
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